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We analyze efficiency of excitation energy transfer in photosynthctic complexes in transient and stationary 
setting. In the transient setting the absorption process is modeled as an individual event resulting in a 
subsequent relaxation dynamics. In the stationary setting the absorption is a continuous stationary process, 
leading to the nonequilibrium steady state. We show that, as far as the efficiency is concerned, both settings 
can be considered to be the same, as they result in almost identical efficiency. We also show that non- 
Mar kovianity has no effect on the resulting efficiency, i.e., corresponding Markovian dynamics results in 
identical efficiency. Even more, if one maps dynamics to appropriate classical rate equations, the same 
efficiency as in quantum case is obtained. 
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I. INTRODUCTION 

Excitation energy transfer in the initial stages of photo- 
synthesis has gained large interest due to coherent beat- 
ings observed in two-dimensional electronic spectroscopy 
experiments on photosynthctic complexe*'^. In light of 
these observations, multiple mechanisms have been pro- 
posed that could lead to improved efficiency of energy 
transfei^^l^^. However, the relevance of the experiments 
(and suggested mechanisms) for the actual processes in 
vivo is still debatecP^lt2il^ as photosynthesis takes place 
in natural conditions of incoherent continuous sunlight 
illumination, while the experiments are conducted by a 
coherent pulsed laser light. 

The process of excitation energy transfer (EET) in- 
volves electronic excitations on pigments and molecular 
vibrations of pigments and nearby protein^^H Proper 
treatment of vibrational degrees of freedom (environ- 
ment) in a description of EET is not trivial, as coupling 
strengths in pigment-protein complexes (PPCs) are such 
that the environmental effects can not be treated per- 
turbatively. While in the limit of weak and strong en- 
vironmental coupling Redfield and Forster theory22] give 
simple and intuitive description of excitation dynamics, 
there are many suggested methods that are trying to 
properly account for environmental effects also in the in- 
termediateregime-'^ Recently, hierarchical equations 
of motiorP^* ^^ * ^ * (HEOM ) gained much popularity in the 
context of EEli^ESlSMil^ as it is formally exact, however, 
at the expense of high numerical efforP^. Also, due to in- 
volved mathematical structure, it offers little insight into 
underlying principles governing the dynamics of EET. 

Two different settings for the EET can be considered. 
In experiments with short laser pulses the excitation 
transfer is just a transient phenomenon - an initial exci- 
tation is either transferred to the target site or dissipated 
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in the environmenl)^. After long time there are no excita- 
tions nor currents present. Such a situation will be called 
a transient setting. In natural conditions though there 
is a constant flux of incoming photons that continuously 
create excitations. After a very short transient time a 
stationary state i s esta blished, the so-called nonequilib- 
rium steady stat^^l^, supporting time-independent en- 
ergy flow. This second situation will be called a station- 
ary setting. 

In this work the focus is on a comparison of a transient 
and stationary setting, in partic ular o n the differences in 
the efficiency of EET. Efficienc)^^!^^! corresponds to the 
probability that the absorption event will result in the 
energy being transported to the target site. We study 
two settings because they are physically relevant, i.e., 
transient case for the pulsed light vs. stationary in the 
case of natural light. Also, because the stationary setting 
is by definition time-independent it enables for an easier 
discussion of the role played by various non-Markovian 
and oscillatory effects. The difference between the ef- 
ficiency in the transient and stationary setting is in all 
relevant situations found to be negligible. Therefore, as 
far as the efficiency goes, the two settings are equiva- 
lent. Not least, it turns out that the stationary setting 
can also have some advantage in terms of computational 
speed over the transient setting where the whole time 
evolution has to be computed. 

We consider various approximations when analyzing 
the efficiency, each providing description at a differ- 
ent level of detail. We start with a generalized quan- 
tum master equation, which provides a complete de- 
scription of EET dynamics, including non-Markovian ef- 
fects due to the interaction with environment. The ker- 
nel for a generalized master equation is obtained from 
the HEOM method. From the generalized quantum 
master equation we obtain the corresponding Markovian 
quantum master equation, and, following the Nakajima- 
Zwanzig formalisnpH, also the corresponding classical 
master equation. We shall show that the efficiency is 
identical in all three cases, i.e., for the HEOM, Marko- 
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vian approximation, as well as for simple classical rate 
equations. Also, main features of the EET dynamics are 
retained at each level of approximation. This result sug- 
gests that simple rate equations might be adequate for 
the description of the processes relevant for the biolog- 
ical function of PPCs provided the calculation of rates 
properly takes into account the underlying quantum me- 
chanics. 



II. MODEL 

Dynamics of excitations in photosynthetic complexes 
can be described at different level of detail, and can be 
either based on derivation from microscopic picture, or 
phenomenological with parameters obtained from experi- 
ments. First we will classify equations of motion (EOMs) 
based on their mathematical structure, ignoring underly- 
ing microscopic model. We will also introduce a formal- 
ism that enables a consistent mapping of EOMs from full 
quantum description to the level of classical rate equa- 
tions. In the following subsection relevant microscopic 
model for PPCs is introduced, providing full quantum 
description of PPCs based on the HEOM method. Note, 
however, that the finding about the equivalence of effi- 
ciencies of EET and the role of non-Markovianity does 
not depend on the specific form of the microscopic model 
used. 



A. Types of EOMs 

Microscopic description of the photosynthetic system 
is given by the total density matrix of the system R{t), 
containing electronic and vibrational degrees of freedom 
(DOF) of pigments and surrounding proteins. Evolution 
of R is governed by the Schrodinger equation 



(1) 



Such complete description is however computationally in- 
tractable due to large number of DOFs. Therefore, the 
total system is usually divided to a relevant (system) and 
an irrelevant (environment) part, with the relevant part 
corresponding to electronic DOFs and the irrelevant to 
the vibrational DOFs. Then the effective EOMs are de- 
rived for the system density matrix p{t) only. The proce- 
dure is formally exact by Nakajima-Zwanzig formalism'^* 
by introducing projection operators for the relevant and 
irrelevant part V and Q that act on the total density op- 
erator, where V is chosen such that p{t) pph — VR{t). 
Projectors satisfy usual relations = V, = Q and 
V ^ Q — H. When the initial state i?(0) and the projec- 
tor Q are such that Qi?(0) ~ 0, the following equation is 
obtained, 



dp(i) 
dt 



= r/c(i 



T)p{T)dT, 



(2) 



which is known as a generalized quantum master equa- 
tion, and contains only the relevant density matrix of 
electronic DOF. However, calculation of the kernel K.{t) 
from the microscopic picture of eq. ([l]) is highly nontriv- 
ial. Nonetheless, for certain cases of system-environment 
interaction, efficient numerical schemes have been devel- 
oped, enabling an exact evolution of system density ma- 
trix. Most frequently used in the context of EET are the 
hierarchical equations of motion (HEOM) , which are also 
used in the present paper. In HEOM the direct evalua- 
tion of memory kernel lC{t) and time-nonlocal evolution 
is circumvented by the introduction of auxiliary opera- 
tors. The details of the method will be given in next 
subsection. 

In certain regimes the time-nonlocal equation ([2]) can 
be simplified by the Markovian approximation in which 
the kernel is taken to be /C = K.5{t), i.e., there are no 
memory effects, resulting in a time-local quantum master 
equation, 



dpjt) 
dt 



= fCpit). 



(3) 



Quantum Markovian eq. ([s]) can also serve as a star- 
ing point for the derivation of the corresponding classical 
dynamics, i.e., equations dictating the evolution of di- 
agonal elements of system's density matrix in a certain 
basis, p = {poo, pn, . . . , Pnn)j that is of populations. The 
corresponding classical generalized master equation is of 
the form 



dpjt) 
dt 



K{t - T)p{T)dT, 



(4) 



and with an additional Markovian approximation a clas- 
sical master equation is obtained. 



dp{t) 
dt 



Kp{t). 



(5) 



Formally, one can derive classical master equation 
from quantum master equation by employing Nakajima- 
Zwanzig formalism, where the projection operators V 
and Q are chosen to project out only dynamics of popu- 
lations (see appendix |A] for details). 

In the present work we shall use the term non- 
Markovian for evolutions governed by a time-dependent 
kernel, eqs. ([2| or Q, while we call evolution Markovian 
if it is determined by a time-local kernel, eqs. ([s]) or ([5|. 



B. Microscopic model 

Here we specify the microscopic model of PPC that is 
usually employed when treating EElEH. The EOMs de- 
rived from the model result in a generalized master equa- 
tion ([2|. We start by separating the Hamiltonian into 



two parts, H = H 



ppc 



Hi, 



where "Hppc corresponds 



to an isolated PPC (electronic and vibrational DOFs), 
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and 'Hint accounts for the electro-magnetic field interac- 
tion (leading to absorption/recombination) and interac- 
tion with other nearby functional units (e.g. reaction 
center). In the following, we will treat dynamics due to 
Hppc exactly, while the effect due to Hint will be treated 
approximately on a phenomenological level. 

Hamiltonian for the isolated PPC is decomposed as 



n 



ppc 



He. 



Hph + H. 



cl — ph ; 



(6) 



with 



N N 

Hci = ^ em|m)(TO| + ^ Vjnn\m){n\, (7) 

N N 

- E = E E f^^ibi^b^^, (8) 

m—l m—1 ^ 

N N 

7^ci-ph - E ^e7-ph = E T^aTi^U + bmi)\fn}{m\, 

m—1 m—1 ^ 

(9) 

where N is the number of pigments in PPC, Hci corre- 
sponds to electronic DOFs within single-excitation man- 
ifold, Hph are phonon DOFs due to pigment and pro- 
tein vibrations, and Hei-ph account for exciton-phonon 
interactions. |m) corresponds to the excitation on the 
mth pigment within the single-excitation subspace, 
is the corresponding on-site energy and Vmn accounts 
for the inter-pigment interaction, b^^^^ and bm£, are cre- 
ation/annihilation operators for the ^th phonon mode 
coupled to the mth pigment, is the frequency of the 
corresponding mode, and the coupling of the excita- 
tion on the mth site to the i^th mode. 

Formal solution of eq. ([T]) for the system density matrix 
in the case of separable initial condition i?(0) = po^Pph^ 
Pint is given by 

p(<) = trpi,,int{ exp[(£ppc -I- Cint)t]pph 8) Pintjpo, (10) 

where Liouvillians C are linear superoperators deter- 
mined by their action on a density matrix, Cp — 
— j^[H,p]. Evaluation of time evolution of p{t) is non- 
trivial already in the case of an isolated PPC as the 
pigment-protein interaction cannot be treated perturba- 
tively. Introduction of interaction Hamiltonian Hint com- 
plicates matters even further, as generally [Hppc, Hint] 7^ 
0. For the isolated PPC, exact nonperturbative method 
has been developed that accounts for the Hci-pii inter- 
action by int roducin g a hierarchy of equations of mo- 
tion (HEOMj^mmi for auxihary DOFs. The HEOM 
method can be considered to be an exact description for 
Lorentzian spectral density and will be used as a starting 
point for various approximations that we explore. Dy- 
namics due to Hint will be taken into account approxi- 
mately by extending resulting HEOMs by effective oper- 
ators obtained from Born-Markov approximation^^*^. 



We assume that each pigment is coupled to an inde- 
pendent phonon bath, where the mth bath has a Drude- 
Lorentz spectral density, Jm(a;) ~ ^ 
ujm^), which is 



(11) 



Spectral density is characterized by a reorganization en- 
ergy Am that specifies strength of the interaction between 
excitons and phonons, and the bath relaxation time 7""'^. 
In high-temperature limit, fcT > /17, which is relevant for 
PPC dynamics at room temperature, HEOMs are of the 
fornpSl 



dpnjt) 
dt 



N 



N 00 



= \^ci-J2 "^0^0 1 ~ E E TT t^j' t^j' 



j=i k=i 



N 



J = l 
N . _ 



(12) 



where pn are auxiliary density matrices, accounting for 
memory effects in evolution, and n is a vector enumerat- 
ing them, n = (77.1,77.2, ...,11^). System density matrix 
corresponds to p(t) = Pn=o{t). Formally, we can rep- 
resent HEOM as linear first order differential equation 
dp{t)/dt — Ap{t), where p contains all auxiliary density 
matrices pn- We will refer to the sparse operator A as the 
HEOM operator. Hierarchy of equations is terminated by 
a criterion J2i ^ -^max, where A^max must be chosen 
such that the memory effects of the evolution are appro- 
priately accounted for. is a shorthand notation for a 
vector differing from n in the jth component, Uj — ^ nj±l. 
Vk = 2TTk/ph are Matsubara frequencies, and complex co- 
efficients Cjk are given by Cjo = Xjjj{cot{l3h'yj /2) — i)/h 
and Cjk = 4Aj7-,j^fe/((j/^ - for ^ > 1, where 

(5 = l/(fcT). We have also introduced a shorthand nota- 
tion Vj = We note that the HEOMs can be repre- 
sented as a generalized quantum master equation ^ with 
the procedure for the memory kernel evaluation given in 
appendix [Aj 

The effect of Hint can be included into HEOM by intro- 
ducing an effective time-local Liouvillian acting on 
system density matrix p{t). The corresponding combined 
dynamics can be obtained by augmenting electronic Li- 
ouvillian Cc\ with the effective interaction Liouvillian, 
£ci -5- Cc\ + i^^, resuhing in a hybr id HEOM-Born- 
Markov set of equations of motiorP^EII. The exact form 
of C^^ that is used for the modeling of absorption (i.e., 
pumping), recombination and transfer of excitation to 
a nearby functional units will be given in the following 
sections. 
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(a) transient case 
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FIG. 1. An example of time evolution of electronic den- 
sity matrix p{t) and of probability rate to the sink js{t), cor- 
responding to the (a) transient, and (b) stationary setting. 
Note that in (b) p{t) approaches a stationary state, having 
time-independent pii{t) and js{t), while in (a) the initial ex- 
citation is either transferred to the target site or lost to the en- 
vironment, resulting in a trivial long-time state. Evolution is 
shown for a dimer system with parameters V ■ 
100cm-\ 7 = 10" " ' 
r = 1 X 10-^ is-^, 
IV I. 



-^is-\ A = 100 cm"', K = 1 X 10-'fs-\ 
Q = 1 X lO"** fs"\ T = 300K (see section 



III. EFFICIENCY 

The efRciency of excitation energy transfer in photo- 
synthesis is the probability that the absorption event 
will result in a transfer of excitation to the target func- 
tional unit, commonly a reaction center. The efficiency of 
smaller functional unit can also be considered, in which 
case it corresponds to the probability that the incoming 
excitation (e.g., due to transfer from an antennae) will 
be transferred to the next functional unit (e.g., a reac- 
tion center). Typical example of such smaller functional 
unit is the Fenna-Mattews-Olson (FMO) complex, which 
acts as a linker between a chromophoric antennae and a 
reaction center. 

Depending on the setting we use, stationary or tran- 
sient, the efficiency has to be defined appropriately. In 
the stationary setting it has to account for the absorp- 
tion (i.e., pumping) and subsequent transfer of excitation 
as a continuous stationary process, while in the transient 
setting one has a time-dependent relaxation dynamics 
from the initial excited state. Stationary setting is suited 
for the description of a PPC under natural light con- 
ditions, where individual absorption events are not re- 
solved, while the transient one can correspond to the case 
of absorption due to short light pulse. In previous studies 



the transient setting has been often employecPHIMIEllEZl 
The following analysis demonstrates that the efficiency in 
the transient and stationary setting is almost identical, 
with small difference only due to effect of absorption on 
the internal dynamics of PPC. 

While the dynamics of p{t) due to phonon bath is ex- 
actly treated by previously introduced HEOMs, we still 
have to specify that will be used to model recom- 
bination of excitation, transfer to reaction center and in 
stationary picture also the absorption (or transfer from 
an antennae). The relevant system state space consists 
of single-excitation space |m) and electronic ground state 
|0). Recombination of the excitation to the ground state 
will be modeled by the operator 

^ / 1 \ 

/:rccomb(p)=2r^ |0)(n|p|7l)(0|--{|n)(n|,p} , 
11=1 ^ ^ 

(13) 

where P is a site-independent recombination rate. Trans- 
fer of excitation to the reaction center is modeled by an 
analogous operator 



Ai„k(p) = 2At |0)(s|p|s 



(14) 



where s denotes a site connected to the reaction center. 
Note that the reaction center is not explicitly included 
in p(t), and £sink only causes transition from sink site to 
the ground state |0). To model absorpti on (or t ransfer 
from antennae), we introduce an operatoi^lSlSlIll] 



'Cabs(p) = 2a \a 



1 



,P} 



(15) 



where a denotes the absorption rate, while a is the site 
that gets excited due to the absorption / transfer pro- 
cess. For simplicity, we have chosen the simplest forms 
of £rocomb, ^sink and £abs, howcvcr, they can be triv- 
ially generalized to linear combination of operators on 
different sites, e.g. transfer to sink from multiple sites or 
absorption on multiple sites. 

We shall now define the efficiency in the two settings, 
stationary and transient. In the transient setting the ex- 
citation is initialized at a special "input" site that we also 
call the absorption site because it is the same site that is 
involved in the absorption process in the stationary set- 
ting, pq = |a)(a|. The system's state is then propagated 
with >Cppc and where 

rcS r- _L r 

-^T '^iccomb I -^sink ■ 



(16) 



Initial excitation decays to the ground state either due to 
recombination or due to transfer to the sink. As there is 
no absorption, the system converges to the ground state 
|0)(0| after long time. Once time evolution of p{t) is ob- 
tained, the transient efficiency can be calculated by in- 
tegrating probability rate of transfer of excitation to the 
sink. 



?7t 



js{t)dt^2K {s\p{t)\s)dt, 



(17) 
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where the probabihty rate of transfer to the sink js fol- 
lows from the expression for >Csink- For an example of 
time evolution in the transient case see Fig. [T^. 

For the stationary setting, the stationary state of the 
system po 
is obtained, where 



p{t — ?> oo) under evolution by £ppc and £g*^ 



In addition, using the final value theorem in a situation 
with a unique nonequilibrium stationary state, the fol- 
lowing useful expression for the stationary state |poo)) 
can be obtained. 



lini zf](z)|po)) 



IPoo))- 



(24) 



S -'-'rccomb i *-'s 



abs 



(18) 



Once we have the stationary state poo the efhciency is 
calculated as a ratio between probability rate of transfer 
to the sink jg and probability rate of absorption event 



3c 



Vs 



ia(oo) 



K {s\ Poo \s) 

a (0| Poo |0) ■ 



(19) 



Time evolution of density matrix populations and proba- 
bility rates for the stationary scenario is shown in Fig. [lb. 
Note that the sole difference between the transient (17) 



For convenience we also introduce a Liouville space no- 
tation, i.e., Hilbert-Schmidt space of operators, with 
\mn)) — \m){n\, {{mn\ = (|mn)))^ and a scalar product 
{{A\B)) — tr{A^B). In the analysis we shall decompose 
/C and the corresponding propagators to various contri- 
butions and observe how each term affects the efhciency 
of EFT. We shall also decompose a time-dependent ker- 
nel to the effective Markovian contribution K, and to the 
non-Markovian contribution IC, 



IC{t) = K.S{t)+K.{t) 



(25) 



and the stationary ( 19 ) setting is in the presence of the 



£abs that causes constant pumping of excitations and the 
appearance of a nonequilibrium stationary state. 

Now we return to the description of dynamics via a 
generalized master equation ([2|, with the kernel lC{t) cor- 
responding to the microscopic dynamics due to £ppc and A. Transient efficiency 
C^^ . The exact form of the kernel is calculated (s ee ap- 
pendix |A| using the HEOM-Born-Markov methocP^Ml. 
Note that derivations require only specific form of C^^ , 
while iZppc can be arbitrary. For such general scenario, we 
show that both efficiencies, transient and stationary, de- 
pend only on the time-integrated kernels /C, while actual 
time-dependence of kernels has no effect on the efficiency. 
Also, the difference between stationary and transient ef- 
ficiency is very small for the typical parameters of FET 
in photosynthesis, so both measures of efficiency can be 
considered equivalent. 

Before analyzing each efficiency measure we introduce 
some common tools that are employed in the analysis. 
For the comparison of efficiencies Laplace transform is 
used, 



where the Markovian contribution corresponds to the 
integrated kernel, IC = K,{t)dt, while the non- 
Markovian contribution is JC{t) — IC{t) — K.5{t). 



In the Laplace picture the transient efficiency is ex- 
pressed as 



?7t —2k lim ((ss I fix (2:) 

z->0 



(26) 



which follows from the properties of the Laplace trans- 
form, and ^It:{z) is the propagator for the kernel /Cx via 
eq. ( 22 ) . Kernel /Cx correspond to the dynamics due to 
pigment-protein interaction C^pc and and effective Liou- 
villians for the transient case, of eq. (16). 



p{z) 



'p{t)dt, 



(20) 



resulting in a Laplace-transformed generalized quantum 
master equation ^ as 



We decompose the kernel /Cx to the Markovian and 
non-Markovian contribution /Cx = /Cx^(0+^T(i)- Prop- 
agators are decomposed accordingly using eq. ( [23| , re- 
sulting in ^t{z) = ^t{z) + nT{z)IC'T{z)nT{z), with 
the obvious notation r2x(2) = ^/{z — /Cx)- Inserting 
this expression into the definition of the transient effi- 
ciency ( 26 ) , we obtain two contributions to the efficiency, 
Vt = Vt~^VTi where 



p{z) ^n{z)p{o) 
n{z) =(z-/c(z))-\ 



(21) 
(22) 



77rp = 2k lim {{ss\il,T^(z)\aa)) 



2^0 



(27) 



where n{z) is the Laplace transform of a propagator, 
and IC{z) is the Laplace transform of a memory kernel. 
Laplace-transformed quantities are indicated by their ar- 
gument z. At several occasions we will need a Laplace 
transform of a propagator resulting from a kernel that 
is a sum of two terms, JC{z) — /Ci(z) + K,2{z). Writing 
^l{z) = l/(z-/Ci(z)-/C2(z)) = r!i(z)/(l-/C2(z)l]i(z)), 
we obtain 



n{z) = ni{z) + ni{z)iC2{z)VL{z). 



(23) 



?7x = 2k Iim((ss|r2x(z)/Cx(z)r2x(2)|aa))- (28) 

2— >o 

The Laplace transform of a non-Markovian kernel van- 
ishes for z = due to lC{t)dt = 0, and can thus be 

approximated for small z as JCt{z) « K,^^ ■ z + C(z^), 
and expression for the non-Markovian contribution as 

77X = 2k lim((ss|nx(z)^^^^zl7x(z)|aa)) (29) 

Identifying the limiting expression limz_j.o zflj^{z)\aa)) as 
the stationary state ( 24 1 and noting that in the absence of 
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absorption it is equal to the trivial ground state, |jOoo)) = 

|00)), as well as ^^^^|00)) = 0, it follows that the non- 
Markovian contribution in the transient case vanishes, 
fiT = 0. The efficiency in the transient case therefore 
depends only on the Markovian kernel /Ct, 



(30) 



i.e., it does not dejjend on the non-Markovianity which 
is all contained in )CT{t). 



B. Stationary efNciency 

For the stationary efficiency the steady state \poc,)) is 
unique and therefore independent of the initial state \po)) . 
Using eq. (24) the stationary efficiency (19) can be writ- 
ten as 



Vs 



lim 



K{{ss\ns{z)\po)) 

h a{mns{z)\po)) 



(31) 



Stationary state |poo)) is also the zero-eigenvector of the 
corresponding Markovian kernel, which is evident by 
observing the stationarity condition for the generalized 
master equation ([2]), 



ICit)\poo))dt^lC\poo)) =0. 



(32) 



Therefore, using eq. (19), we can equivalently write the 



efficiency with the Markovian propagator only. 



?/s = lim 



K{{ss\ns{z)\po)) 



Oa((00|Os(z)|po)) 



(33) 



Similarly as in the transient case, the efficiency does not 
depend on the non-Markovian part ICs (t) . To obtain the 
stationary efficiency we therefore only need /Cs and not 
the full kernel ICs{t) = ICsS{t) + JCs{t). We are going to 
write JCs as a sum of a transient Markovian kernel ICt, 
the absorption Liouvillian £abs, and the rest. 



abs 



ph 



(34) 



where (as we shall see small) term IC^^^^ arises due to the 
non-commutativity [£abs, 'Cci-ph] 7^0. Expression for the 
77s can now be further simplified using eq. ( 23 ) , by writing 



the Markovian propagator for the stationary case as a 
sum of propagators for the transient case and the rest, 

fls{z) = nT{z) + ^T{z){£ahs+ICl^s)^s{z), where we also 
used the fact that the Laplace transform of a Markovian 
kernel, being a delta function in time, is equal to the 
kernel itself. Inserting this expression into the numerator 
of eq. (33 1, we obtain 



Vs 



lim 



n z{{ss\nT{z){C^hs + IC^^,s)^s{z)\po)) 



z((00|17s(z)|po)) 



(35) 



where we have taken into account that the stationary 
state for the transient setting is trivial, i.e., only ground 
state is occupied, \miz^Qllss\z^T{z)\pQ)) = ((ss|00)) = 0. 
After inserting the >Cabs from eq. (15), the expression for 



the stationary efficiency becomes a sum of the transient 
efficiency and a correction due to the absorption. 



?7s = »?T + ?/A- 



(36) 



The correction due to the absorption can be expressed as 



T^A = 2k lim 



AT 
m,n— 1 



ph. 



((ss|rjT(z)/Cabsl"in))((mn|rjg {z)\aa)) 



ph 



with the propagator fig (z) = (z — /Ct — /Cabs) 



(37) 
The 



details of the calculation can be found in appendix |Bj 

Numerical calculations in the following sections show 
that the difference between the stationary efficiency and 
the transient efficiency t^a is very small for the micro- 
scopic model of PPC considered. Therefore, for practical 
applications, they can be considered to be the same. Ob- 
serve that if one starts with a Markovian description of 
PPC dynamics, eq. ([3|, adding absorption term £abs to 
obtain a stationary setting, the efficiency difference jyA is 
exactly zero. 

In the analysis above, we have considered transient and 
stationary efficiency in the case of dynamics described by 
a generalized master equation ([2]). As the efficiency only 
depends on the corresponding Markovian kernel /C, the 
dynamics under time-local quantum master equation ([s]) 
results in the identical efficiency of EET, 77s = ?7g and 
?7t — Vt- That is, a detailed time-dependence, e.g., non- 
Markovian oscillations in density matrix elements, has no 
direct effect on the efficiency. 

For classical master equation ([4| and ^ an analogous 
analysis can be conducted, where the effective rates must 
be calculated from the corresponding effective Liouvil- 
lians for the transient or stationary case of eq. ( 16 ) or 
(18 1. Thus, similarly as in quantum case, for classical 



EOMs the efficiency also depends only on the classical 
Markovian kernel K. Moreover, if one projects a quan- 
tum master equation to a generalized classical master 
equation in an appropriate basis, and so that the effec- 
tive Liouvillians translate to the corresponding effective 
rates, the efficiency of EET is the same in both cases. 
Therefore, mapping of dynamics from non-Markovian 
quantum description to Markovian classical description, 
IC{t) — > /C — > K{t) — > K, does not affect the efficiency, 
i.e., all four efficiencies are exactly the same. Mapping of 
time-independent Markovian quantum master equation 
with kernel /C to the corresponding time-dependent non- 
Markovian classical master equation with kernel K{t) 
is done by employing the Nakajima-Zwanzig formalism, 
where a projection operator is chosen such that the di- 
agonal elements of system density matrix represent the 
relevant subsystem. See appendix [A| for details. 

The above findings about the equivalence of the two 
efficiency measures and on the absence of non-Markovian 
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FIG. 2. Time evolution for the transient setting in a dimer 
system and different levels of approximation: exact HEOM 
method (pii(t), solid line; full kernel ICT{t)), Markovian ap- 
proximation (Pii(t), dashed line; ICt) and classical Markovian 
rate equations dot-dashed line; Kt)- Evolution of the 

population of the input site pii(t) (initial state is po = |lXl!) 
and of the matrix element of a non-singular part of the mem- 
ory kernel /Ci2,i2(i) (real and imaginary part) as well as of the 
classical non-Markovian kernel Ki2{t) is shown. Two differ- 
ent reorganization energies are used, in (a) A = 10cm~^ and 
in (b) A = 100 cm~^. Other parameters are V = 100 cm~^, 
e = 100cm"\ 7 = 10"^fs"\ T = 300K. Note that even thou 
time dependence is different, the efficiency rjT is the same in 
all three cases. 



effects do not rely on the specific form of the microscopic 
model, e.g., on the exact form of the spectral density 
J{uj). The value of the efficiency itself of course does 
depen d on m icroscopic parametera^ltlil^like the spectral 
densitji2I122Hll]^ however, rjs and ?7t are affected in exactly 
the same way. 

In the following, we shall calculate the transient and 
stationary efficiencies and the corresponding populations 
dynamics numerically using the HEOM formalism for the 
specific microscopic model described in previous section. 
Equivalence of efficiency measures will be demonstrated 
as well as mapping of a generalized master equation to a 
corresponding Markovian classical master equation using 
the Nakajima-Zwanzig formalism. 



IV. DIMER SYSTEM 

As an illustrative example we analyze the case of a two- 
site PPC. The interaction with environmental phonons is 
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FIG. 3. Calculated efficiency of excitation transfer for the 
dimer system at different values of reorganization energy A. 
Transient efficiency t^t is shown (solid line, left axis) and the 
absorption contribution in the stationary case j^a (dashed 
line, right axis). Observe that the difference ?7a is of order 
~ 10~^ and therefore, for practical purposes, 77T ~ '7s- Pa- 
rameters are the same as in Fig. [2] 



treated exactly within the HEOM formalism from which 
the time-dependent kernel /C(i) is explicitly evaluated 
(see appendix |A] for more details) . Transient efficiencies 
77t and stationary efficiencies X]^ are calculated for a range 
of reorganization energies A, demonstrating that the dif- 
ference between the efficiency measures 77A is negligible. 
For the transient case, mapping of dynamics from non- 
Markovian generalized master equation ([2| to the cor- 
responding classical master equation ([5| via Nakajima- 
Zwanzig formalism is also demonstrated. 

The total Hamiltonian for the dimer is of the form ([6]) 
for two sites, N = 2. The relevant parameters of the elec- 
tronic Hamiltonian "Hci are the on-site energy difference 
e = £2— El and the inter-site interaction strength V = Vi2. 
Each site is coupled to an independent bath with identi- 
cal parameters A„ — A and 7„ — 7. First site is an input, 
a — 1, while the second is connected to the reaction cen- 
ter, s = 2. We consider parameters within ranges typical 
for PPCs. Thus, if not stated otherwise, the tempera- 
ture is T = 300 K, bath relaxation time 7 = lO^'^fs^^, 
while recombination rate F = 2 x 10^^ fs~^ and sink rate 
K — 2 X 10"'* fs""'^. For the stationary case the absorp- 
tion rate is taken the same as the relaxation rate, a = F, 
however, the efficiency 77s (and the corresponding differ- 
ence from the transient case 77A) does not depend on the 
actual choice of a. 

In Fig. [2] exact time evolution of the input site popu- 
lation pii(i) (starting from the initial state po — |1)(1|) 
is shown for two values of reorganization energy A in the 
transient setting. Emergence of incoherent dynamics is 
evident as the reorganization energy is increased, result- 
ing in a faster decay of coherent oscillations seen at short 
times. Time dependence of matrix element /Ci2,i2(i) of a 
non-singular part of the kernel ( A3 1 is also shown, where 
we use a short notation /Ci2,i2(0 = [I^ns{t)]i2,i2- Note 
that other non-zero matrix elements of IC{t) also decay 
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on a comparable time scale. Markovian dynamics for the 
integrated memory kernel K, is also calculated. Compar- 
ing the time dependence of population on the site 1 for 
the exact non-Markovian evolution and the Markovian 
approximation we can see that the Markovian evolution 
results in a faster decay of coherent oscillations. In addi- 
tion, oscillations for the non-Markovian and Markovian 
case, although different, are such that the efficiency is 
the same in all three cases. 

We also mapped dynamics from quantum Markovian 
master equation with kernel K, to the corresponding 
classical generalized master equation using Nakajima- 
Zwanzig formalism, obtaining a time-dependent classi- 
cal kernel K{t). Time dependent rate K2i{t) is shown 
in Fig. [2] The dynamics of classical populations p{t), 
eq. Q, is identical as the dynamics of the diagonal el- 
ements of for quantum Markovian case because the 
mapping is exaclP^J. Integrating time-dependent kernel 
K(t), classical time- independent master equation with 
kernel K and population dynamics p{t) is obtained, 
eq. ([5). Again, the decay of initial oscillatory dynamics 
is evident when approximating dynamics of p(t) with the 
corresponding classical Markovian dynamics p{t). Here 
we emphasize that while the populations dynamics un- 
der mapping p{t) — )■ p{t) are all different, the resulting 
efficiency of EET is identical. 

For the stationary case ([Ts]) exact steady state is ob- 
tained by finding the zero-eigenvector of a sparse HEOM 
operator A from which the corresponding stationary ef- 
ficiency ryg is calculated. We have calculated rj^ for a 
range of reorganization energies A. The results are shown 
in Fig. [Sj For the parameters taken, the contribution due 
to the absorption is of order of ^ 10~^, indicating that for 
the purpose of analysis of PPC efficiency, both settings 
can be considered to be equivalent. 
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FIG. 4. (a) Time evolution of populations pmm{t) for the 
transient evolution of the FMO model and the initial state 
Po = |1X1|- Dashed lines represent corresponding Markovian 
evolution In the inset, a short-time dynamics for pn, 

P22, Markovian dynamics Pj^j^ and P22 (dashed black line) and 
classical Markovian dynamics Pj , P2 (dash-dotted red line) is 
shown, (b) Real part of selected elements of a non-singular 
part of the time-dependent kernel K.{t) for the FMO model. 
Parameters of simulation are 7"^ = 166 fs, A = 35cm~^, 
K = 2 X 10-*fs-\ r = 2 X 10-'^fs-\ r = 300K. 



V. FENNA-MATTEWS-OLSON COMPLEX 

The Fenna-Matthews-Olson complex (FMO) is often 
considered in the studies of PPC as its structure is well 
known and a large number of variou s experimental stud- 
ies has been conducted on itPEI^. We have used the 
HEOM method to calculate the transient and stationary 
efficiencies and to evaluate the memory kernel. 

We considered a 7-site FMO model, with electronic 
Hamiltonian "Hei as specified in Ref. 031 Site 1 is con- 
sidered as an input site, while the site 3 is connected 
to the reaction center. Each site is interacting with an 
identical phonon bath with parameters A = 35cm ~^, 
= 166 fs, already used in previous studie^^I^Hii], i^e- 
combination and sink rates were taken from Ref. I32|, with 
r = 2 X 10-6 fs-i and K = 2 X lO-^fs"!. 

Time evolution of populations p{t) for the transient 
case is shown in Fig. |4^; dashed lines in addition show 
Markovian dynamics ]}{t) due to K,. For Markovian dy- 
namics the initial oscillatory behavior of populations is 
not as pronounced as in a dimer, while at latter times 



Markovian and non-Markovian evolutions result in al- 
most the same site populations. In the inset a short time 
dynamics for sites 1 and 2 is shown. Additionally, popula- 
tions p{t) for classical master equation, resulting from the 
mapping K. K{t) K, are also shown (dash-dotted 
line) . At latter times populations p{t) approach the val- 
ues corresponding to the full quantum master equation 
p{t). Note again that the efficiency does not depend on 
the level of approximation. 

To provide some insight on the time-dependent kernel 
that is relevant in a PPC, few entries of the non-singular 
part of kernel Knsit) are shown in Fig.[4]D. Plotted entries 
contribute to the decay of coherences between site 1 and 
other sites. Note that the element /Ci2,i2 has a larger 
decay time than other elements, which can be related to 
the fact that the site 1 is most strongly coupled to the 
site 2 in 'Hei- 

For the used parameters the difference between the sta- 
tionary and the transient efficiency is less than the esti- 
mated error of the calculation due to the truncation of 
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the HEOM hierarchy at A^max = 8, with the efficiencies 
estimated to be rys « t^t ~ 0.962, and the truncation 
error being of the order ^ 10^"*. We have estimated the 
truncation error by observing the convergence of station- 
ary efhciency for calculations with HEOM truncation at 
level up-to A^max = 12. For practical purposes the effi- 
ciencies ryx a-nd rjg can thus be considered to be the same. 

We note that due to the equivalence of efficiency mea- 
sures in stationary and transient setting one might choose 
to calculate the one that is easier to evaluate. Stationary 
efficiency turns out to be more convenient in that respect, 
i.e., its calculation can be faster than calculating whole 
time evolution of Pn{t), as it requires only the calculation 
of a stationary state poo, being the eigenvector with the 
zero eigenvalue of the HEOM operator A. Taking advan- 
tage of the sparsity of A and using iterative eigenvalue 
algorithms we obtain the stationary state. 



the mapping from full non-Markovian quantum picture 
to the corresponding classical picture is properly treated. 
This is consisten t with previous essentially classical de- 
scriptions of EE'S^S'^. The procedure for obtaining clas- 
sical rates from a microscopic model is however highly 
non-triviapsMi ^i^q coupling to environmental DOFs 
can not be treated perturbatively. Starting from the ex- 
act quantum description, we use Nakajima-Zwanzig for- 
malism and the Markovian approximation to obtain a de- 
scription on a relevant subspace. Thus, we obtain, as far 
as the efficiency is concerned, equivalent descriptions of 
dynamics at different levels of detail. This also provides 
a straightforward way of comparing different approxima- 
tions. 



Appendix A: Evaluation of memory kernel from HEOM 



VI. CONCLUSION 

We have studied two physically relevant settings of 
EET, namely, a transient situation in which the ini- 
tial excitations decay to zero, and a stationary set- 
ting in which constant pumping causes the system to 
converge to a nonequilibrium stationary state. Com- 
paring efficiency between the transient and stationary 
setting we observed that the difference is very small 
for the exact description, while it is exactly zero for 
Markovian equations, like, e.g., the Lindblad equa- 
tion. Equivalence of efficiency measures in transient 
and stationary case validates findings about efficiencies 
in PPCs based on o bserving transient dynamics from 
specific initial state^^HniMllSHlTl Therefore, the mech- 
anisms leading to higher efficiencies established in the 
transient p icture, suc h as environment-as sisted quantum 
transpor 

^mmmmm and supertransfeiESIlSl, are also 
relevant in the case of incoherent light illumination, com- 
plying with the findings of Ref. 16, 

In both settings, transient and stationary, the ef- 
ficiency does not depend on a (time-dependent) non- 
Markovian part of the kernel. Same result also holds 
for the description with the classical rate equations. Ad- 
ditionally, if one obtains classical rate equations from full 
quantum description by employing Nakajima-Zwanzig 
formalism, the resulting efficiencies are the same as in 
full quantum description. The only feature of the dy- 
namics that is not reproduced by the classical or quan- 
tum Markovian equations is the initial oscillatory motion 
of populations being present for pure initial states. While 
physical relevance of evolutions from pure initial quan- 
tum states for th e in vivo process of photosynthesis is 
questionabl^i^EII, this oscillatory motion, being present 
or not, does in no way affect the corresponding efficiency 
of EET, i.e. an approximate dynamics without oscilla- 
tory character results in identical efficiencies. 

This suggests that the rate equations are adequate for 
the description of EET in biological processes, as long as 



We shall evaluate the memory kernel for system den- 
sity matrix based on the hierarchical equations of motion. 
We are treating HEOM as evolution dp{t)/dt = Ap{t), 
where .4 is a time-independent linear sparse operator de- 
fined by eq. (12 1. Memory kernel is obtained from A 



by projecting out auxiliary degrees of freedom ^,15.^0 us- 
ing Nakajima-Zwanzig formalisni'^^, resulting in a sum of 
singular and non-singular contribution 

JC{t)^K,sm+^ns{t), (Al) 

which can be evaluated in Schrodinger picture as 



JCs =VAV, 
ICnsit) ^VAg{t)QAV. 



(A2) 
(A3) 



V and Q are projectors to pn=a and Pn^Oi respectively, 
and g{t) is the propagator for the irrelevant part Qp, 



g{t)=exp{QAt), 
being a solution of a differential equation, 

dg{t) 



dt 

with the initial condition 



QAGit), 



5(0) = /. 



(A4) 



(A5) 



(A6) 



Numerically, one can evaluate each column of g{t) in- 
dividually as d[g{t)]j/dt = QA[g{t)]j, where [g]j is the 



jth column of Q. However, direct evaluation of eq. (A5) 



becomes intractable when the number of sites A^ and 
the hierarchy truncation level A'max is increased, even 
if sparsity of operator A is taken into account. The 
number of auxiliary matrices in HEOM is given bj0 
A^tot = (A^ + Afmax)!/(A^Wmax!)i whilc cach auxiliary ma- 
trix has A^^ elements. Dimensionality of HEOM operator 
A is thus A'^ = N'^Ntof Obtaining Q{t) requires numer- 
ical solution of Nj^^ differential equations for vectors with 
A'_4 components. 
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A complete solution of Q (t) is however not required to 
obtain the kernel. This is evident if we introduce the 
operator 



nit) = Ag{t)QAP, 

which is calculated by integrating 
dTZ{t) 



(A7) 



dt 



= AQAg{t)QAV = AQn{t), (A8) 



with the initial condition 

7^(0) = AQAV. 



(A9) 



From the operator TZ, a non-singular part of the kernel 
is obtained as 



iCnsit) = vnit). 



(AlO) 



The main advantage of calculating TZ{t) instead of Q{t) 
is that due to a projection in eq. (AlO), only columns 



of TZ(t) that correspond to relevant degrees of freedom 
have to be calculated when obtaining the kernel ICns{t), 
i.e., only solution of N"^ differential equations is required 
instead of N'^Ntot- 

When obtaining classical generalized master equa- 
tion Q from quantum master equation mapping of 
a time-independent kernel /C to a time-dependent kernel 
K(t) is also obtained by the procedure introduced above. 
Kernel /C assumes the role of A, and operators V and Q 
are projectors on the diagonal and off-diagonal elements 
of density matrix, respectively. 



the time-dependent picture, where it corresponds to the 
time-integral of the the sink site population for the evo- 
lution with fix and initial state |00)). The last term is a 
contribution due to the effect of phonon environment on 
the absorption. With the corresponding stationary state 
we can write 



1 lim ((gg|^T(^)/C|^t|poo)) 



a z- 



Poo)) 



(B3) 



Alternatively, we shall rewrite the expression for 
eq. ( B3 ) , so that it does not explicitly dependent on the 
stationary state |poo)) and the absorption rate a. 

First, we note that K^Xs does not contribute to the 
absorption, i.e. 



ph 

abs 



= 0. 



(B4) 



This can be seen if one writes expression for singular and 
non-singular part of the kernel from eqns. ( A2 ) and (A3 ) 



for all DOFs (system and environment), i.e., full Liouvil- 
lian C = £ppc -I- takes the role of A, and projector 
is a projector acting on a full density matrix, resulting 
in VR{t) = p{t) (g) pph. Acting with the corresponding 
expressions for the kernel on the state |0G)) and taking 
into account that only >Cabs acts on the ground state, we 
obtain 



/C,,s|00)) =P£P|00)) = £abs 
/C„.,s(i)|00)) ^VCg{t)QU^J> 



0. 



(B5) 
(B6) 



Appendix B: Contributions to stationary efficiency 

In this appendix, we start with the expression for the 
stationary efficiency in Laplace picture, eq. (351, and 
show that it can be written as a sum of transient effi- 
ciency TjT and contribution due to the effects of phonon 
environment on the absorption. We start be rewriting 



absorption Liouvillian from eq. (15) as 



/labs -2a(|aa))((G0|- |00))((00|) 



N 



2a ^(|0to))((0to| -t- |toO))((Oto|). 



(Bl) 



m— 1 



Inserting this expression into eq. (35), we obtain 



r]s —2Klmi{{ss\QT^{z)\aa)) + 2k \im{{ss\QT{z 

2:— >-0 2— >0 

— lim 

((00|17s(z)|po» 



(B2) 



We identify the first term as the efficiency for the tran- 
sient case r/rp, eq. (27), while the second term is zero. 



This is evident if we convert the limiting expression to 



The absorption effect of a non-singular part is zero due 
to QV — 0, while the absorption effect of a singular part 
corresponds to the absorption Liouvillian >Cabs- The ab- 
sence of the absorption term in /Cabs then follows from 
the comparison with the decomposed stationary kernel 
in eq. ([34]). 

Next, we start by inserting identity into 



the numerator of eq. ( B3 ) , obtaining 



TjA = K lim ((ss|0t(2:)/C 



ph 



N 



abs / J 
m,n— 1 



((mn Poo)) 

Poo)) 



Note that the resulting sum does not contain the ground 
state term |00)) due to (B4). The stationary state jpoo)) is 
expressed with the stationary propagator f2s , which is de- 



composed according to the eq. ( |23| ) as ^s{z) ^ jtg 



rig {z)Cahs^s{z) with r2g"(z) — (z — /Ct 
serting decomposed propagator fig (z) into the numerator 
of eq. ( |B7[ ), we obtain 



■ph. 



ntiz) + 



In- 



{{mn\poo)) 

Poo)) 



{{mn\zV.f' {z)£i,bs^siz)\po)) 
z^s{z)\P())) 



— lim 



(B8) 



=2 lim ((mn|r2g (z)\aa)). 

z— >o 
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In the first line, we have taken into account that 
the stationary state of fig (z) is a ground state, i.e. 

lim^_j.o((TOn|zf^S i^)\Po)) = {{mn\00)) = 0. In the second 
hnc, we have inserted the expression for C^hs from eq. 



(Bll, taking into account vanishing limiting expression 
limz^Q{{mn\n^\00)) = and cancel common terms in 



the numerator and the denominator. Plugging eq. (B8) 



back into (B7|, we obtain the efficiency in Laplace pic- 
ture, 



N 

riA = 2k lim > 

m,n— 1 



{{ss\nT{z)IClljmn)){{mn\nf\z)\aa)). 

(B9) 

Note that while the expression for tja does not explicitly 
depend on the absorption rate a, /Cabs could in principle 
depend on it. We have however verified numerically that 
/Ci^bs does not depend on a for the microscopic model of 

PPC used in this work by calculating /C^bs fo'" ranging 
over several orders of magnitude. 



iH. Lee, Y.-C. Cheng, and G. R. Fleming, S cience 316, 1462| 
13 (2007) 

^G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. Mancal, 
Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming, |Nature| 
Q 446, 782 (2007) 
^E. CoUini, C. Y. Wong, K. E. Wilk, P. M. G. Curmi, P. Brumer, 

and G. D. Scholes, Nature 463, 644 (2010) 
*G. Panitchayangkoon, D. Hayes, K. A. Fransted, J. R. Caram, 
E. Harel, J. Wen, R. E. Blankenship, and G. S. Engel, |Proc.| 
I Natl. Acad. Sci. U. S. A. , 1 ( 2010) 

°P. R ebentrost, M. Mohseni, and A. Aspuru-Guzik, J. Phys!] 
I Che m. B 113, 9942 (2009) 
°A. Ishiz aki, T. R. Calhoun, G. S. Schlau-Cohen, and G. R. 
Fleming, iPhys. Chem. Chem. Phys. 12, 7319 (2010)1 



''L. A. Pachon ancTP. Brumer, J. Phys. Chem. Lett. 2, 2728 (2011_ 
*A. Ishizaki and G. R. Fleming, Proc. Natl. Acad. Sci. U. S."A 

r 106, 17255 (2009) 

^M. B. Plenio and S. F. Huelga, New J. Phys. 10, 113019 (2008) 
^"P. Rebentrost, M. Mohseni, I. Kassal, S. Lloyd, and A. Aspuru- 
Guzik, New J. Phys. 11, 033003 (2009) 
'^-'^M. Mohseni, P. Rebentrost, S. Lloyd, and A. Aspuru-Guzik, J. 
PChem. Phys. 129, 1 74106 (2008) 

^■^F. Carus o,1^ W. Chin, A. Datta, S. F. Huelga, and M. B. 

Plenio, J. Che m. Phys. 131, 105106 (2009)) 
13 J. Wu, F. Liu, Y. Shen, J. Cao, and R. J. Silbey, |New J. Phys.| 
r 12, 105012 (2010) 



L. A . Pachon and P. Brumer, |Phys. Chem. Chem. Phys. 14^ 
I 10094 (2012) 

"P. Brumeran d M. Shapiro, |Proc. Natl. Acad. Sci. U. S. A. 109,| 
I 19575 (20"l2)1 



1®I. Kassal, J. Yuen-Zhou, and S. Rahimi-Keshari, |j. Phys. ChemT] 

Lett. 4, 362 (2013) 
^'T. Mancal and L. Valkunas, New J. Phys. 12, 065044 (2010) 
'^^Y.-C . Cheng and G. R. Fleming, An nu. Rev. Phys. Chem. 60] 
\ 241 ( 2009) 

^«F. FassioU, A. Olaya-Castro, and G. D. Scholes, [x Phys. Che"im] 

Lett. , 3136 (2012) 
^"M. Tiersch, S. Pope scu, and H. J. Briegel, [Phil. Trans. R. Soc.| 

A 370, 377 1 (2012)] 
■=^L. A. Pachon and P. Brumer, 'Phys. Rev. A 87, 022106 (2013)" 
^^V. May and O. Kiihn, Charge and Energy Transfer Dynamics in 

Molecular Systems, 3rd ed. (Wiley- VCH, 2011). 

23a. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 23411l| 

(2009) 

24p. Huo and D. F. Coker,| j. Chem. Phys. 133, 184108 (2010)| 
25G. Ritschel, J. Roden, W. T. Strunz, and A. Eisfeld, |New J.| 

Phys. 13, 113034 (2011). 
'^°D. P. S . McCutcheon and A. Nazir, |j. Chem. Phys. 135, 11450f] 
> (2011)1 

^'P. Nalbach, D. Braun, and M. Thorwart, |Phys. Rev. E 84,| 
041926 (2011) 

^» Y. Tanimura and R. Kubo, |J. Phys. Soc. Jpn. 58, 101 (1989)| 
29Q. Shi, L. Chen, G. Nan, R.-X. Xu, and Y. Yan, |j. Chem. Phys!] 

130, 084105 (2009) 
^"J. Zh u, S. Kais, P. R ebentrost, and A. Aspuru-Guzik, | J. Phys.| 

Chem. B 115, 1531 (2011) 
^^A. G. Dijkstra and Y. Tanimura, New J. Phys. 14, 073027 (2012)^ 
3^0. Kreisbeck, T. Kramer, M. Rodriguez, and B. Hein, J. Chem, 
[Theory Comput. 7, 2166 (2011), 

Ritz, S. Park, and K. Schulten,|j. Phys. Chem. B 105, 82 59| 

(2001) 

•'^H.-P. Breuer and F. Petruccione, The theory of open quantum 

systems (Oxford University Press, USA, 2002). 
35a. W. Chin, A. Datta, F. Caruso, S. F. Huelga, and M. B. 

Plenio, New J. Phys. 12, 065002 (2010) 
36a. Shabam^M. Mohseni, H. Rabitz, and S. Lloyd, |Phys. Rev. E| 
p86, 1 (2012)1 



I 



■^'S. Jesenko a nd M. Znidaric, >Jew J. Phys. 14, 093017 (2012) 
38d. Manzano, [PLoS ONE 8!'e57 041 (2013)| 



39c. Kreisbeck and T. Kramer,rj7Phys. Chem. Lett. 3, 2828 (2012)] 
4°A. KoUi, E. J. O'Reilly, G. D. Scholes, and A. Olaya-Castro, [j] 
r~Chem. Phys. 137, 174109 (2012) 

^M. del Rey, A. W. Chin, S. F. Huelga, and M. B. Plenio, [j] 
[~~Phys. Che m. Lett. 4, 903 (2013) _ 
^Mapping of quantum master equation with kernel K to classi- 
cal generalized master equation with kernel K{t) via Nakajima- 
Zwanzig formalism is exact only when the initial state po is diag- 
onal, i.e., no coherences are present. Otherwise, inhomogeneous 
terms have to be included in classical master equation. 
*3j. Adolphs and T. Renger, Biophys. J. 91, 2778 (2006) 
^'^T. Brixner, J. Stenger, H. M. Vaswani, M. Cho, R. E. Blanken- 
ship, and G. R. Fleming, Nature 434, 625 (2005) 
■^^S. Lloyd and M. Mohseni, New X Phys. 12, 075020 (2010)] 
46e. N. Zimanyi and R. J. Silbey7 p: Chem. Phys. r33, 144107| 
(2010) 

*'J. Briggs and A. Eisfeld, Phys. Rev. E 83, 4 (2011) 

Cao and R. J. Silbey, J. Phys. Chem. A 113, 13825 (2009)| 
'^^J. Wu, F. Liu, J. M a, R. J. Silbey, and J. Cao, |J. Chem. P'^^ 
I 137, 174111 (2012)1 



